The Effect of Disorder in an Orbitally Ordered Jahn- Teller Insulator 



Sanjeev Kumar^, Arno P. Kampf^, and Pinaki Majumdar^ 
^ Institute of Physics, Theoretical Physics III, Center for Electronic Correlations and Magnetism, 
University of Augsburg, D-86135 Augsburg, Cermany 
^ Harish- Chandra Research Institute, Chhatnag Road, Jhusi, Allahabad 211 019, India 

We study a two dimensional, two-band double-exchange model for Cg electrons coupled to Jahn- 
Teller distortions in the presence of quenched disorder using a recently developed Monte-Carlo 
technique. In the absence of disorder the half-filled system at low temperatures is an orbitally 
ordered ferromagnetic insulator with a staggered pattern of Jahn- Teller distortions. We examine the 
finite temperature transition to the orbitally disordered phase and uncover a qualitative difference 
between the intermediate and strongly coupled systems, including a thermally driven insulator to 
metal crossover in the former case. Long range orbital order is suppressed in the presence of disorder 
and the system displays a tendency towards metastable states consisting of orbitally disordered 
stripe-like structures enclosing orbitally ordered domains. 

PACS numbers: 71.10.-w, 75.47.Lx, 81.16.Rf 

February 4, 2008 



I. INTRODUCTION 

The undopcd pcrovskite manganites, for example 
LaMnOa, arc orbitally ordered antiferromagnetic insula- 
tors at low temperature. The magnetic order is of A-type 
with ferromagnetic planes coupled antiferromagnetically 
in the transverse direction. The orbital order is of C- 
type, anti-correlated with the magnetic order 0, and is 
accompanied by long range ordering of the Jahn- Teller 
(JT) distortions of the MnOg octahcdra. Upon hole dop- 
ing, e.g. in Lai_2;Ca^Mn03, the system evolves via a fer- 
romagnetic insulating phase towards an orbitally disor- 
dered (OD) ferromagnetic metal, which exhibits colossal 
magnetoresistance (CMR) near the Curie temperature 
Q . Quenched disorder has been widely acknowledged as 
a crucial ingredient for the understanding of hole-doped 
manganites • It is the source of a variety of phe- 
nomena, both in real materials as well as in model cal- 
culatioiis, like the coexistence of metallic and insulating 
phases [3, laB, [1,3, thermally driven metal to insulator 
transitions |l0llllLll2l[T3l |. and the suppression of charge 
order in half-doped systems 0, 0| . 

Although the 'fame' of the manganites rests on the 
CMR observed near optimal hole doping (x ^ 0.3) the 
low doping materials also display a rich variety of phe- 
nomena 0, 01 . Hole doping leads to progressive loss 
of orbital order (00), a reduction of the transport gap, 
increasing isotropy in the effective magnetic exchange, 
and an insulator (I) to metal (M) transition 0, ITsL . 
Beyond transport and thermodynamic indicators, NMR 
and neutron scattering experiments have revealed an in- 
homogeneous - possibly phase separated - state in the 
low doping regime [23, llj. The interplay of doping, 
disorder, and thermal fluctuations thus presents an in- 
triguing problem even in the insulators. In typical hole 
doped manganite compounds such as REi_j,AEa;Mn03, 
where RE (AE) denotes a trivalent (divalent) rare-earth 
(alkaline-earth) ion, the disorder arises from the differ- 
ence in the ionic radii of the RE and AE ions and hence. 



is tied to a variation in the hole concentration, vanishing 
for the stoichiometric compositions at a; = and x — \. 
Conceptually, however, it is useful to disentangle the ef- 
fects of increasing hole density and increasing disorder on 
the 00-OD and IM transitions. In that spirit, this pa- 
per focuses on the effect of disorder and varying electron- 
lattice coupling at a; = in a two-orbital model for the 
Mn Cg-electrons coupled to JT lattice distortions and to 
the t2g-derived core spins. The combined effect of disor- 
der and hole doping has been briefly reported elsewhere 
0. A model system, which includes disorder remain- 
ing at a; = is experimentally relevant for systems with 
isoelectronic substitutions such as REi-^REj^MnOa [2^ . 

The manganites are complex materials involving the 
the interplay of spin, orbital, lattice, and charge degrees 
of freedom. In the cubic perovskite structure, the 5-fold 
degeneracy of the Mn-d levels is lifted by crystal fields 
leading to two manifolds of three t2g and two eg orbitals. 
The energetically low-lying t2g levels contain three elec- 
trons and give rise to a well localized S = 3/2 spin. The 
Cg electrons are itinerant and interact with the distor- 
tions of the MnOg octahedra via the JT coupling and 
with the core spins via Hund's rule coupling jl^. The 
full complexity of these interactions is hard to handle 
in a realistic two or three dimensional situation. Let us 
quickly survey the attempts to understand the x — 
state before describing our approach. 

Different points of view have been reported regard- 
ing the relative importance of the inter- and intra- or- 
bital Hubbard interactions and the JT interaction in the 
manganites [23L 01 ■ Although the electron-electron (e- 
e) interactions are large in these systems, it has been 
suggested that the parent insulator can not be consid- 
ered canonical Mott insulator Surprisingly, 
the orbitally ordered insulating groundstate at a; = 
can be understood from a variety of starting points. 
For strong Hubbard interactions a spin-orbital t-J model 
has been derived and studied within mean-field approx- 
imation [26l [2^ leading to an orbitally ordered ground- 
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state. The orbitally ordered groundstate also emerges 
from Hartree-Fock band-structure calculations consider- 
ing the Mn-3d and 0-2p orbitals [H i^. Monte-Carlo 
studies of models including JT coupling and ignoring the 
local Coulomb repulsion also explain the OO-I state and 
tend to describe the system as a JT insulator ■ In fact 
e-e and JT interactions are known to have qualitatively 
similar consequences and instead of competing they re- 
enforce each other [sj ■ More recently it was shown that 
a Fermi surface nesting instability at weak coupling can 
also be the source for an orbitally order ed g roundstate 
in the presence of Hubbard interactions jsj. It is also 
known that a more realistic modelling of the lattice tak- 
ing into acount the individual oxygen displacements of 
MnOg octahedra and their cooperative effects, can by it- 
self lead to a staggered JT distorted, and hence orbitally 
ordered, groundstate [s^. Retaining all the complica- 
tions of the real materials is a difficult task, particularly 
when we wish to advance beyond mean field theory. The 
principal purpose of the present study instead, is to clar- 
ify the influence of thermal fluctuations and disorder on 
the orbitally ordered insulating groundstate. Here we ig- 
nore the e-e interactions and the cooperative character of 
the lattice distortions. The similar effects of e-e and JT 
interactions in determining the nature of the groundstate 
suggests that the effect of disorder may also be similar 
for the two cases. 

The remainder of the paper is organized as follows. In 
Section II we specify the model and briefly describe our 
Monte-Carlo simulation technique. The results are pre- 
sented in Sec. Ill, starting with the non-interacting sys- 
tem and subsequently introducing the JT coupling and 
quenched disorder. The non disordered zero temperature 
limit is independently analyzed by using variational cal- 
culations, while the clean strong coupling limit studied 
via an effective classical model. 



II. MODEL AND METHOD 

We consider a two-band model for itinerant eg elec- 
trons on a square lattice. The electrons are coupled to 
JT lattice distortions, t2g derived S = 3/2 core spins and 
quenched disorder as described by the Hamiltonian: 

{l])cr i 

i i i 

Here, c and are annihilation and creation operators 
for Cg electrons and a, (3 are summed over the two Mn- 
Cg orbitals d^2_y2 and c?3;j2_^2, which are labelled (a) 

and (6) in what follows, denote the hopping ampli- 
tudes between Cg orbitals on nearest-neighbor sites and 
have the cubic pcrovskite specific form: t^^ = t^^ = t, 



ti, = tl, ^ t/3, ti, = ^ ti, =jj„ t/V3, 

where x and y mark the spatial directions |30j . Disorder 
is modelled by random on-site potentials with equally 
probable values ±A, which couples to the local electronic 
density n.i . The Cg-electron spin is locally coupled to the 
t2g spin Si via the Hund's rule coupling Jh- The elec- 
tronic spin is given by erf = J2aa' 4aa^acr'^'iaa' , where 
arc the Pauli matrices. A denotes the strength of the 
JT coupling between the distortion Qi — {Qix,Qiz) and 

the orbital pseudospin rf = Y."^ 4aa^a0^tP'^- K is a 
measure of the lattice stiffness, and ^ is the chemical 
potential. We set t = 1 = K as our reference energy 
scale. The JT distortions and the <2g derived core spins 
are treated as classical variables, and we set |S| = 1. The 
present study is restricted to half-filling and we explore 
the variation in the parameters A and A in addition to 
the temperature T. 

For further simplification we adopt the limit Jh /t ^ 
oo, which is justified and frequently used in the context 
of manganites 0, ^3- In this limit the electronic spin 
at site i is tied to the orientation of the core spin S;. 
Transforming the fermionic operators to this local spin 
reference frame leads to the following 'spinless' model for 
the Cg electrons: 

ap 

(ij) i 

-A^Q,.T. + |5:Qf. (2) 

i i 

The new hopping amplitudes have an additional depen- 
dence on the core spin configurations and are given by: 

— = cos — cos-;f +sm — sm-;^ e . (3) 

tafj 2 2 2 2 

Here, 9i and 0, denote polar and azimuthal angles for the 
spin Si . From now on the operator Cia {c\^ ) is associated 
with annihilating (creating) an electron at site i, in the 
orbital a with spin parallel to S^. 

The model given by Eq. (2) is bilinear in the electronic 
operators and does not encounter the problem of expo- 
nentially growing Hilbert space, since all many-particle 
states can be constructed from Slater determinants of the 
single-particle states. The difficulty, however, arises from 
the large phase space in the classical variables Q and S. 
At zero temperature the problem reduces to finding the 
spin and lattice configurations {Si,Qi} that minimize 
the total energy. The energies for a limited number of 
periodic structures in Q and S can be compared analyti- 
cally. This method is not assured, though, to lead to the 
true groundstate and often requires additional physics in- 
sight for identifying periodic structures, which are prime 
candidates for the groundstate. The finite temperature 
properties are not accessible in this manner, since they 
necessarily require the electronic energies and wavcfunc- 
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tions in non-periodic structures for the Q and S vari- 
ables. Further complications arise in the presence of dis- 
order and even the groundstate may not belong to the 
subspace of periodic configurations of the spin and lat- 
tice variables. Two controlled methods are available at 
finite T: (a) dynamical mean-field theory (DMFT) cor- 
rectly captures the strong coupling physics 0, |33 but 
is unable to handle spatially correlated inhomogeneous 
states, which may emerge in the presence of quenched 
disorder, while, (&) exact diagonalization based Monte 
Carlo (ED-MC) calculations [l^ can handle both ther- 
mal fluctuations and disorder, but they are severely size 
limited. 

Here, we will use a scheme which is closely related to 
ED-MC and employ the travelling cluster approximation 
(TCA) in order to overcome the small size limita- 
tions. This approximation is based on the observation, 
that the effect of a local update on spin and lattice vari- 
ables does not 'propagate' long distances via the elec- 
trons, and the energy change involved in such a move 
can be computed by constructing a cluster Hamiltonian 
around the reference site. This drastically reduces the 
computational cost and allows access to lattice sizes of 
~ 1000 sites. For the electronic properties in the ther- 
mally equilibrated system we use ED for the full system. 
The details of this computational scheme have been dis- 
cussed previously [s^. While TCA forms the backbone 
of the present study and is used to map out the phase di- 
agram of our model Hamiltonian in the parameter space 
of A, A, and T with detailed real-space information, we 
also analyze the limiting cases using variational calcula- 
tions, or an effective classical model for the JT lattice 
distortions. Apart from providing a comparison to the 
TCA results, this also enables us to establish a more 
transparent physical picture of the numerical results. 

III. RESULTS 

A. Orbital order in the groundstate 

We begin by describing the simplest limit of the Hamil- 
tonian in Eq. (2), which is the non- interacting electron 
system in the absence of disorder. Since the magnetism 
is purely double-exchange driven, leading to a ferromag- 
netic (FM) groundstate, a fully polarized core spin state 
is used for the calculation of the electronic dispersion 
relation. For A = the spectrum is straightforwardly 
obtained by Fourier transformation to momentum space: 

^»-^E ''^''■"4«- > (4) 

k 

This gives 

H = J2 ^"Mk) d^du^ (5) 

k,a/3 

with eQ/3(k) = —2t^pCos{kx) — 2i'^pCos{ky). Diagonaliz- 
ing this Hamiltonian we obtain 
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FIG. 1; (Color online) (a) Density of states for the non- 
interacting Hamiltonian. Solid (dashed) line is for the 
(E^) band and the total density of states is shown as the 
dotted line, (b) Filling of the d^2_y2 and d^^-2_^-2 orbitals as 
a function of the chemical potential. The solid line is the total 
filling 71. (c) Fermi surfaces for the non-interacting model at 
n = 1. The solid (dashed) segment of the Fermi surface arises 
from the band E^ (E^)- The three nesting wave- vectors are 
indicated by the arrows. 



where e='=(k) = eaa(k) ± £66 (k). 

Fig. 1(a) shows the band- resolved and the total density 
of states (DOS) corresponding to the dispersion given in 
Eq. (6). The electronic DOS is given by: 

n 

where En denotes the eigenvalues of the Hamiltonian. 
For the band-resolved DOS we use or E^^ in place of 
Fig. 1(b) plots the filling of the d^2_y2 and (1^22 -r'^ 
orbitals as a function of the chemical potential. These are 
simply the groundstate expectation values of the number 
operators dlj^^ and 46^^k6- 

While the occupancy of the two orbitals is equal at 
n ~ 1, d.j.2_y2 is preferred for < ?i < 1. Fig. 1(c) shows 
the Fermi surface at electron filling n = 1, which corre- 
sponds to /i = 0. From panel (a) it is clear that both 
bands arc partially filled at /i = and the Fermi surface 
consists of contributions from both bands. While each of 
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the bands provides a segment of the Fermi surface, 
nested with wavevectors (7r,0) and (0, tt), there is also 
interband nesting with wavevector (tt, tt). The existence 
of the nesting wavevectors and a therefore a divergent 
non-interacting susceptibihty induces an ordering insta- 
biUty in the presence of interactions. For example, the 
inclusion of an infinitesimal electron-lattice coupling. A, 
is expected to lead to a lowering in the electronic energy 
by stabilizing the lattice patterns, which are compatible 
with a selected nesting wavevector. Similar arguments 
have been invoked earlier for the same band dispersion 
to show that the presence of weak e-e interactions leads 
to orbital ordering at half filling • 

In order to confirm this in the context of electron- 
lattice interactions, we perform a variational calculation 
to search for the lowest energy state within a restricted 
set of possible groundstates. As discussed earlier, the 
T = problem amounts to finding the spin and lattice 
configurations {Si, Q^}, which minimize the total energy. 
Since the groundstate is known to be ferromagnetic, we 
are left with the problem of determining the minimum- 
energy lattice configuration {Q} only, for which we set 
up a restricted variational calculation. We consider vari- 
ational states of the type Qi = Q e"! ""' and compare 
energies for q = (0,0), (0, tt), (tt, 0) and (tt, tt). The 
variational parameters are the magnitude of the distor- 
tion Q and the orientation of Q in the — Qz plane 
parameterized by an angle C, with tan((^) = Qx/Qz- For 
this restricted set of lattice configurations the Hamilto- 
nian matrix is reduced toa2x2or4x4 form, which 
is diagonalized to evaluate the total energy and to con- 
struct an approximate T = phase diagram. The matrix 
for the uniform distortions (q = 0) is 



(a) 0I 



2 

Q = 







(b) 



1 - 



1^2 




FIG. 2: (Color online) (a) Schematic picture of a staggered 
distortion pattern for tlie JT distortion vectors Q. The ar- 
rows indicate the directions of Q, with an arbitrary angle ^. 
The groundstate corresponds to i^a ~ 7i"/2, and (^b = 3tv/2. 
(b) The variation of the magnitude of the distortion shown 
in (a) with A. The circles are obtained from the variational 
calculation and the squares are from the Monte-Carlo sim- 
ulations, (c) Density of states for A = 0, 1 and 1.4 in the 
variational groundstate. Lattice sizes of up to 2000 x 2000 
sites are used in the variational study. 



eii(k)-AQz ei2(k)-A(5j; 
ei2(k)-AQ,. e22{k) + XQ, 



and the matrix for q = (0, 7r),(7r, 0) or (tt, tt) is 



/eii(k) ei2(k) 

£12 (k) £22 (k) 

—XQz —XQx 



"XQz —XQx 
—XQx XQz 

eii{k + q) ei2(k + q) 

ei2(k-t-q) e22(k + q) 



The results of the variational calculation arc summa- 
rized in Fig. 2. A staggered pattern of JT distortions 
with an arbitrary orientation angle C is schematically 
shown in Fig. 2(a). The minimum energy is obtained 
for (a = 7i'/2 and = 3tt/2. The origin of purely Qx- 
typc distortion patterns is tied to the specific structure 
of the hopping parameters tap and especially to the sign 



difference between tj^^ and t^f^. 



For tip = 



f^p the min- 



imum energy configuration corresponds to a staggered 
distortion pattern with ~ ^[2tab/itaa ~ tab)] and 
Cs = tt + Ca- For t^ij — t'^i^ = t/-\/3 this leads to (^a = 7r/3 
and (b = 47r/3, and for t^f^ = = —t/^/3 the corre- 
sponding angles are 27r/3 and Stt/B. The relative sign 
difference between and t^^^ in cubic pcrovskitcs, thus 
leads to a conflict regarding the orientation angles for 



the directions of Q on the two sublattices and (a = 7i'/2, 
= 37r/2 result as a compromise. The Qa:-type distor- 
tion patterns can also be motivated from the structure 
of the effective classical model at strong JT coupling, 
which is again related to the perovskite specific hopping 
parameters. Since the lattice distortions are coupled to 
the orbital pseudospin, Q^^-type ordering in the JT dis- 
tortions is accompanied by Ta;-type staggered order in the 
orbital sector. 

Fig. 2(b) shows the the magnitude Qmin of the dis- 
tortion corresponding to the minimum energy state as 
a function of the electron-lattice coupling. For weak A 
there is a 'BCS' like instability leading to an exponen- 
tially small Qmin and a correspondingly small gap in the 
DOS [33. The window 0.0 < A < 1 roughly corresponds 
to this weak coupling regime and displays a slow rise in 
Qmin with A. For 1 < A < 2, the intermediate cou- 
pling regime, Qmin increases rapidly crossing over to the 
strong coupling asymptote, Qmin oc A, for A > 2. The 
squares mark the data points obtained from TCA, which 
are discussed later. Fig. 2(c) shows the DOS for the lat- 
tice structure of panel (a) in the weak and intermediate 
coupling regime. 

At strong coupling it is reasonable to start from the 
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atomic limit and to assume that the electrons are site- 
localized by strong lattice distortions [s^. In that case 
the total energy per site, which includes the electronic 
and the elastic energy, is given by 

i?tot =-A|Q| + y|Q|2. 

Minimizing the total energy to obtain the optimum dis- 
tortion Qmin, one finds Qmm — X/K. This is precisely 
the observed behavior for A > 2. Such a self-trapped ob- 
ject is usually referred to as a polaron and the associated 
energy is termed single-polaron energy Ep = /2K. In 
this paper, a 'polaron' always refers to a 'static polaron' 
since the lattice is treated in the adiabatic limit. The 
qualitative differences between intermediate and strong 
coupling will become apparent in the finite temperature 
studies using TCA. 
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FIG. 3: (a) Temperature dependence of the lattice structure 
factor Dq ((tt, tt)). (b) Variation of the orbital-ordering tran- 
sition temperatures Too with A. The squares are the TCA 
results and the solid line is a guide to eye. The strong coupling 
asymptotic form is also shown by the dashed line. 



B. Monte Carlo results 

1. TCA in the clean limit 

The finite T problem is solved using the TCA method 
described earlier. In this case an unrestricted search is 
performed with respect to the spin and lattice configu- 
rations and the temperature is reduced step by step to 
obtain the groundstate. Fig. 3(a) shows the tempera- 
ture dependence of the q ~ (tt, tt) = qo component of 
the lattice structure factor, 

y 

Here and below {.■■)th denotes the average over thermal 
equilibrium configurations. All the TCA results are ob- 
tained on a 24 X 24 lattice and a 4 x 4 travelling cluster. 
Dgiq^o) is a measure of the staggered ordering tendency 
of the JT distortions. Since the lattice distortions are 
coupled to the orbital pseudospin the same tendency is 
transferred to the analogously defined orbital structure 
factor Dr{qo). Therefore -Dq (qo) serves as an indicator 
for the staggered ordering of both, the lattice distortions 
and the orbital pseudospin. The point of inflection in 
the temperature dependence of Dq (qo) is taken as the 
orbital ordering transition temperature Too- Fig- 3(b) 
shows the variation of Too with A. The symbols are TCA 
results and the solid line is a guide to eye. The extrapola- 
tion of the solid line to A < 0.8 is based on the existence 
of a BCS-like instability at weak A due to Fermi-surface 
nesting as already encountered in the variational calcula- 
tions. In the weak-coupling regime Too is proportional 
to the energy gap Aq in the DOS. The dashed line in- 
dicates the result of a strong coupling expansion, which 
leads to the proportionality Too oc t'^/Ep. The details 
of the strong coupling expansion are discussed in the ap- 
pendix. The non-monotonic behavior of Too is thus nat- 
urally understood by merging the weak and strong cou- 
pling limits. The experimental results on REMnOa for 



the accessible range of decrease in mean ionic radii show 
an increase in Too H^- A reduction in ionic radii leads 
to a decrease in the bandwidth, which translates to an 
effective increase in X/t in our model. This suggests that 
the relevant regime for the strength of JT coupling should 
be A < 1.8. A more direct estimate for the value of A is 
obtained by comparing the ratio Ac /Too between our 
results and the experiments. For LaMnOa, Aq ^ 0.5eV 
[13 and Too 700X leading to Aq/Too ^ 8. This 
value is approximately reproduced for A ~ 1.2. In real 
materials additional interactions may lead to a renormal- 
ization of the effective JT coupling jsj, therefore these 
numbers should be considered only as a rough estimate. 

Since the orbital order is tied to the ordering of the JT 
distortions, it is useful to track the temperature evolution 
of the distribution function for the lattice distortions 



P{Q) = 




where | | denotes the magnitude of the local lattice dis- 
tortions for an equilibrium configuration. The tempera- 
ture dependence of P{Q) is shown in Figs. 4(a) and 4(b) 
for moderate and strong JT coupling strengths, respec- 
tively. The low T distributions are peaked at Q = Qmin, 
suggesting an unimodular structure at T = and sup- 
porting the variational results for the groundstate. The 
symbols represent the TCA data and the solid lines follow 
from simple arguments for the expected structure of the 
distribution function P{Q): The single-site energy for an 
electron self-trapped by a distortion of magnitude Q is 
given hy E = -XQ + K/2 Q^. The probability that this 
site has a distortion with magnitude Q at a given tem- 
perature T is then given by P{Q) oc exp[ — {E — Eo)/T], 
where Eq is the groundstate energy correspnding to Q = 
Qmin- This leads to P{Q) cx exp[ - [Q - Qminf I'^T], 
assuming that the functional form of the energy does 
not change for finite T. Using the Qmin data of the 



6 




P(Q) 



0.3 



N(a)) 
0.1 





r (c) i 






4^ 


'1 f 


1 



4 
0) 




FIG. 4: (Color online) Temperature evolution of the distri- 
bution function for the magnitude of the lattice distortions 
P{Q) for moderate (a), and strong (b) electron-lattice cou- 
pling strengths. Symbols are the Monte-Carlo data and the 
solid lines correspond to the naive expectation based on the 
atomic limit. The density of states for the parameters corre- 
sponding to (a) and (b) is shown in (c) and (d), respectively. 



for system sizes ranging from N = 8'^ to N = 32^, a de- 
tailed finite size scaling has not been performed. The ob- 
served qualitative difference between the moderate and 
strong coupling systems raise the possibility, that the or- 
bital order to disorder transitions in the two cases may be 
quaUtatively different. Experimentally it is known that 
the transition in LaMnOa is first order in nature and be- 
comes second order for PrMnO.s and NdMnOa , which are 
smaller bandwidth materials |23, l40l l4ll| . Although the 
origin for a first order nature of the transition in LaMnOs 
is suggested to be anharmonic coupling between JT dis- 
tortions and volume strain [4^ , the bandwidth variation 
may also be playing a role. Figs. 4(c) and 4(d) show the 
electronic DOS for the same parameter values as used 
in (a) and (b), respectively. The low T DOS is gapped, 
consistent with our variational results. The gap fills up 
for T ~ Too for weaker A but persists to much larger T 
for strong coupling. The origin of the gap in the DOS is 
therefore very different for the two cases and finds sup- 
ports in the structure in P{Q). Indeed, if the distortions 
arc well formed even at large T, as is the case for strong 
couphng, the DOS is gapped due to 'self-trapping'. On 
the other hand a gap in the DOS for moderate coupling 
is tied to the existence of orbital order, which in turn 
requires a long range ordered pattern for the lattice dis- 
tortions. Note that the sizes of both Qmin and the gap 
in the DOS obtained in the TCA calculation match very 
well with the variational calculations (see Fig. 2(b)). 



TCA results these functions are plotted as solid lines 
in Figs. 4(a) and 4(b). The TCA results show that, 
P{Q = 0) is finite for T > Too at A = 1.4; thus some 
of the electrons can apparently escape from JT-distorted 
cages upon heating. This suggests the possibility of an 
insulator to metal crossover and will be discussed in the 
next section. The naive analysis described above does not 
apply for weak coupling, as is apparent from the large de- 
viations of the firm lines from the symbols , infact, even 
the width at low temperature of the distribution function 
is not captured by the corresponding Gaussian function. 
For A = 2 however, the electrons are well trapped by the 
self- generated JT distortions even at T ~ STqo- The 
Gaussian functions with width oc Vt are reasonable fits 
to the TCA data, implying that A = 2 is close to the 
regime where single-site analysis is valid. This clarifies a 
crucial difference between the weak and the strong cou- 
pling systems despite the fact that the orbital order in 
the groundstate is the same in both limits. In the strong 
coupling regime the distortions exist even at high tem- 
peratures and the orbital ordering transition corresponds 
to the alignment of the distortion vectors at T ~ Too- 
This is analogous to the ordering in a spin model with 
magnetic moments of fixed magnitude. The mechanism 
for orbital ordering at weak to moderate coupling relies 
on a simultaneous generation and ordering of the lattice 
distortions. 

Although the stability of these results has been checked 



2. The effect of disorder 

Although the presence of disorder is usually tied to 
the hole-doped materials, it can also be realized in com- 
positions like REi^jjRE'yMnOa. Moreover it is useful to 
study the effect of disorder without involving the compli- 
cations of finite hole density. Therefore, given the pres- 
ence of orbital order in the half-filled clean system, we 
now include the effects of quenched disorder. Although 
the Monte-Carlo results provide us with the full finite-T 
information, we first present our results on the disorder 
effects at low T. In this section we focus on the thermo- 
dynamic quantities, which are averaged over ~ 10 realiza- 
tions of quenched disorder. Fig. 5(a) shows the staggered 
component of the lattice structure factor at low T as a 
function of the disorder strength A. The critical value 
Ac for the disappearence of orbital order in the ground- 
state is estimated from these data and further confirmed 
by the T dependence of DQ{qo). Z?q (qo) ~ 0(1) is con- 
sidered to indicate orbital order and £'Q(qo) ^ 0(1/7V^) 
implies an orbitally disordered phase. A non-monotonic 
dependence of Ac on the JT coupling strength is ob- 
served as shown in Fig. 5(b). The similarity between 
Ac (A) and TooW (see Fig. 3(b)) suggests that the ef- 
fects of quenched disorder and thermal fluctuations are 
similar in weakening the long range orbital order. 

For further details on the disorder driven orbital order 
to disorder crossover, we show the disorder evolution of 
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FIG. 5: (a) Lattice structure factor Dq{{-k,-k)) at T = 0.01 
as a function of tiie disorder strength A. (b) Critical dis- 
order strength Ac required to spoil the orbital order in the 
groundstate as a function of A. 



the distribution function P{Q) for the JT distortions for 
two representative values of A in Fig. 6(a) and 6(b). 
These results are very similar to the T dependence of the 
distribution function (see Figs. 4(a)-(b)). Interestingly 
disorder acts as a delocalizing agent for the weak coupling 
system, leading to a fraction of sites with very weak JT 
distortions, hence delocalizing a fraction of electrons from 
the self generated JT traps. The effect of the disorder on 
P{Q) is barely visible at large A and P{Q) is not affected 
in crossing over from the orbitally ordered to the orbitally 
disordered state. Figs. 6(c) and 6(d) show the effect of 
disorder on the DOS for the same values of parameters as 
used in panels (a) and (b). Our earlier suggestion that 
the gap in the DOS is related to the orbital ordering 
for weak electron-lattice coupling is confirmed. For A = 
0.8 a pseudogap feature replaces the clean gap in the 
DOS. The strong coupling DOS is only slightly affected 
by disorder and in particular the gap survives even in the 
orbitally disordered phase. This is again similar to the 
temperature dependence of the DOS for the weak and 
strong coupling limits. 

The effect of disorder is also different for the finite T 
properties between weak and strong coupling. Fig. 7 
shows the temperature dependence of £>Q(qo) for vari- 
ous A at moderate (panel (a)) and strong (panel (b)) JT 
couplings. The T = value of £'Q(qo) reduces with in- 
creasing disorder in both cases. The detailed real-space 
analysis of this reduction will be presented later, but 
it is worth mentioning already here that orbitally dis- 
ordered stripe-like structures emerge in the disordered 
system at low temperatures, specially in the strong cou- 
pling regime. While the inclusion of disorder affects the 
saturation values for DQ(qo) in a similar manner for both 
A = 1.4 and A = 2.0, Too scales show an interesting con- 
trast between the two coupling strengths. Too reduces 
gradually for A = 1.4 with increasing disorder strength 
and eventually both the saturation value of Dq{c{o) and 
Too approach zero as the long range orbital order in the 
groundstate is lost. On the other hand. Too does not 




P(Q) 
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FIG. 6: (Color online) Disorder dependence of the distribu- 
tion function P{Q) for moderate (a), and strong (b) electron- 
lattice coupling strengths. The density of states for the pa- 
rameters corresponding to (a) and (b) is shown in (c) and (d), 
respectively. 



change with disorder for A = 2.0 and drops abruptly to 
zero close to Ac, when the system crosses over to the 
orbitally disordered phase. 

For further confirmation on the difference between the 
weak and strong coupling systems we compute the optical 
conductivity g(a;) using the Kubo formula with the ex- 
act eigenstates |43| . The resistivity p is approximated by 
the inverse of a(ujmin), where tOmin = 20t/N ^ 0.03t is 
the lowest reliable frequency scale for 17(0;) calculations 
on our TV = 24^ system. Fig. 7(c) shows the resistiv- 
ity for A = 1.4 as a function of A for different values of 
disorder strength. A sharp upturn in the resistivity is 
observed upon cooling, which is clearly due to the onset 
of orbital ordering. For all values of A < Ac a change of 
slope in the resistivity is observed, indicating a thermally 
driven I-M crossover. For A > Ac, the system retains 
dp/ dT > for all T, indicating a disorder-induced metal- 
lization of the insulator. While the Anderson localization 
effects may not allow a metallic state in lower dimensions, 
the I-M transition induced by quenched disorder via the 
destruction of orbital order is expected in three dimen- 
sional systems too. The observation of a disorder-induced 
metallic state is not new and has been reported before 
in the context of charge ordering 0, ^| . The tempera- 
ture dependence of the resistivity for A = 2 is shown in 
Fig. 7(d). This strong coupling system displays a more 
robust insulating behavior. The effect of thermal fluctu- 
ations and/or disorder is barely visible on the resistivity. 
Contrary to the weak-coupling system, the onset of or- 





0.1 0.2 0.3 0.2 0.4 
T T 

FIG. 7: Temperature dependence of the lattice structure 
factor Dq{{h, tt)) for varying disorder strengths at (a) A = 1.4 
and (b) A — 2.0. Resistivity p in units of tij-Ke} , as a function 
of temperature for (c) A = 1.4 and (d) A = 2.0 for the same 
values of A as in (a) and (b). 



bital ordering is not reflected in p{T). 

We observe a monotonic decrease of Too scales upon 
increasing disorder in the experimentally relevant regime 
of moderate JT coupling. In real systems however, it is 
difficult to isolate the individual effect of quenched dis- 
order. Isoelectronic rarc-carth substitution series such as 
in Lai_j,NdyMn03, involve a simultaneous change of ef- 
fective bandwidth and disorder p3 | . The increase in Too 
upon increasing y in this series is understood mainly from 
the increase in the average ionic radius but changes in y 
are necessarily accompanied by variations in disorder. 

One of the advantages of the Monte-Carlo method is 
the access to detailed real-space information. This is 
especially useful for determining the mechanism lead- 
ing to the loss of orbital order upon increasing disorder 
strength. In Fig. 8 we show the local spatial correla- 
tions of the JT distortions, as represented by the quan- 
tity Cq = J Qi ■ Qi+5, where i5 is summed over the 
nearest neighbors of site i. Therefore Cq ~ —Qmin 
represents a region of staggered orbital ordering while 
Cq ~ corresponds to an OD neighborhood. Panels 
(a)-(c) show the groundstate pattern for Cq at moder- 
ate strength of the JT coupling for three typical disor- 
der realizations. In some cases orbitally disordered fil- 
amentary structures are observed. These OD domain 
walls separate the orbitally ordered regions, and a tt- 
phase shift appears in crossing from one 00 region to 
the other. The patterns for Cq at strong coupling are 
shown in panels (d)-(f). Again the data arc presented 



FIG. 8: (Color online) Local spatial correlations Cq (see 
text) of the Jahn- Teller distortions at T = 0.01 for three typ- 
ical realizations of quenched disorder. Panels (a)-(c) corre- 
spond to A = 1.6, A = 0.5, and panels (d)-(f) are for A = 2.0, 
A = 0.6. 

for three typical realizations of quenched disorder. The 
tendency to form filamentary structures of orbitally dis- 
ordered regions is even stronger for strong JT coupling. 
A careful analysis of these structures shows that these 
are not 'unique' groundstates of the system and should 
be understood as metastable states. Even for a fixed re- 
alization of disorder a different starting configuration of 
the lattice variables in the Monte-Carlo simulations leads 
to a different filamentary structure. Moreover, in most 
cases these states are higher in energy as compared to the 
states with no domain walls. Upon doping the 00 insu- 
lator, similar stripe-like structures of OD regions appear 
as true groundstates and serve as the prefered locations 
for the doped holes 0. While the earlier study involved 
a combined effect of quenched disorder and hole-doping, 
here we propose that disorder itself leads to the existence 
of metastable low energy configurations with filamentary 
OD regions enclosing 00 domains. The disorder in our 
model couples to the total charge density and in order 
to understand how the effect is transferred to the orbital 
degrees of freedom, we analyze the effect of single impu- 
rity in the orbitally ordered groundstate. Considering a 
repulsive impurity, the impurity site leads to a slightly 
lower value of charge density. Since the band dispersion 
is such that ^ Ub for n 7^ 1, the sites with n < 1 has a 
finite value of . These finite values upset the purely 
T^-type order. 



IV. CONCLUSIONS 

In summary, we have studied a two-band double- 
exchange model at half filling for Cg electrons coupled to 
Jahn- Teller distortions in the presence of quenched dis- 
order in two dimensions. The existence of long range or- 
bital order and the accompanying staggered order in the 
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Jahn- Teller distortions in the groundstate is verified by 
means of variational calculations, a strong coupling anal- 
ysis, and real-space Monte-Carlo simulations. The quali- 
tative difference between the moderate and strong Jahn- 
Teller coupling systems is reflected in the nature of the or- 
bital ordering transitions in the two regimes. A thermally 
driven crossover from an insulator to a metal is observed 
for weak to moderate coupling, in contrast to a robust 
insulating phase at strong coupling. Quenched disorder 
destablizes the orbitally ordered groundstate and leads to 
the appearance of metastabe states with orbitally disor- 
dered domain walls separating the orbitally ordered re- 
gions. These metastable, inhomogeneous states in the 
insulating half-filled system are important for the under- 
standing of structures and phase transitions in the hole 
doped systems 0- 
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APPENDIX: EFFECTIVE CLASSICAL 
MODEL 

In the limit of large JT coupling, the electrons are es- 
sentially site-localized and the atomic limit (tap = 0) is 
a good starting point. The electronic eigenvalues for a 
single-site problem are given by ±A|Q|. Only the lower 
electronic level is occupied for the half-filled case. Mini- 
mizing the sum of electronic and lattice energies for the 
optimum value of the magnitude of the JT distortion Q 
leads to Qmin = X/K. While the magnitude of the JT 
distortion is determined locally by a balance between the 
elastic and the electronic energy, there is no prefered di- 
rection for the orientation of Q. This degeneracy is lifted 
by including a finite hopping amplitude for electrons and 
leads to an effective classical model for the orientations 
of the JT lattice vectors Q. The derivation is straight- 
forward and the main steps are outlined below: 

Consider a two site problem with H' = Hq + V 



-A J2 iQ- 



i=l,2 

V = ^{tafj c[aC2l3 + h.C.) 
a/3 



.^f} + flQ.P (7) 



(8) 



The eigenvalues of Hq are ±A|Qi| with i = 1,2. The 
eigenvectors \ip±) are site-localized, and consist of linear 
combinations of the two orbitals at a site. The eigenvec- 
tor for the eigenvalue — A|Qi| (A|Q.j|) corresponds to r 
pointing parallel (antiparallel) to Q, which has the full 




FIG. 9: Lattice structure factor Dq ((vr, vr)) as a function of 
temperature for the effective strong-coupling model compared 
with the TCA results at A = 3.0. Dq{{tt,-k)) is normalized to 
its saturation value for the TCA data. 



rotational symmetry of the plane. Explicitly, the eigen- 
vectors are: 



COSI- 



|0), (9) 



4 lo),(io) 



where d = tan^^ (Qix/Qiz) denotes the orientation of 
the JT distortion vector Qi in the Qx-Qz plane. 
The total energy of this 2-site system is 



K. 



IQ 



1^ 



(11) 



which is minimized by |Qi| = X/K independent of i. 
The two electron groundstate of Hq thus corresponds to 
having one electron at each of the two sites with a lattice 
distortion |Q| = X/K, total energy Eq = and the 
groundstate wavefunction \ipo) = \tp-)i (8) \tp-)2- 

Turning on the hopping part of the Hamiltonian V, 
the pcrturbative correction to the energy is given by 



ae=j: 



(V'm|ta/3(CiQC2/3 + /l.C.)|V'o) 



En — E„ 



(12) 



where, IV'm) denotes an excited state of Hq in the two- 
electron subspace. Explicit evaluation leads to 



Ai? = Jx sin ^1 sin ^2 + Jz cos Ci cos C2 , 

+Jm (sin Ci cos C2 + cos Ci sin C2) , (13) 

with coupling constants Jx = ^{taahb + ^at) ' 



Jz 



j^y — 2 — 



tlb) ' and J„ 



[tab{taa ^ tbb)] ■ 
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This energy acts as an effective classical Hamiltonian for 
the ordering of the distortion vectors Q. 

For a comparison with the results obtained from the 
TCA, we show the temperature dependence of £>Q(qo) 
for the effective model and TCA in Fig. 9. These re- 
sults are obtained on the same lattice sizes as used for 
the TCA calculations, i.e. 24 x 24. The TCA curve is 
normalized by its T = value. The coefficients used in 
Heff have an additional factor of 1/2 because the sys- 
tem is still paramagnetic when the orbital ordering takes 
place and the hopping amplitude t has to be replaced 



by t\/(l + Si.Sj)/2, which becomes l/-\/2 in the para- 
magnetic phase. The Too scales obtained from H^ff 
match very well with the TCA scales for A > 3 and begin 
to deviate significantly only for A < 2 (see Fig. 3(b)). 
The purely Q^-type staggered ordering in the JT distor- 
tion vectors is also observed within the effective classical 
Hamiltonian. This can be easily understood by looking 
at the values of coupling constants Jx, Jz and Jm- Again 
the difference in sign between the values for and J^, 
which originates from the sign difference between t^j, and 
^ab' turns out to be crucial. 
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